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ABSTRACT 

Optical stellar intensity interferometry with air Cherenkov telescope arrays, composed 
of nearly 100 telescopes, will provide means to measure fundamental stellar parameters 
and also open the possibility of mo del- independent imaging. In addition to sensitivity 
issues, a main limitation of image recovery in intensity interferometry is the loss of 
phase of the complex degree of coherence during the measurement process. Neverthe- 
less, several model-independent phase reconstruction techniques have been developed. 
Here we implement a Cauchy-Riemann based algorithm to recover images from simu- 
lated data. For bright stars (m v ~ 6) and exposure times of a few hours, we find that 
scale features such as diameters, oblateness and overall shapes are reconstructed with 
uncertainties of a few percent. More complex images are also well reconstructed with 
high degrees of correlation with the pristine image. Results are further improved by 
using a forward algorithm. 

Key words: techniques: high angular resolution - instrumentation: interferometers 
- stars: imaging - stars: fundamental parameters 



1 INTRODUCTION 

Even though Stellar Intensity Interferometry (SII) was 
abandoned in the 1970's, there has been a recent interest in 
reviving this technique, mainly due to the unprecedented 
■ (u, v) plane coverage that future imagin g air Cherenkov 
telescope (IACT) array s will provide (|CTA Consortium! 
l20ld : iMaier et ail 120091 ). The possibility of probing stars 
at the sub-milliarcsecond scale and visible wavelengths 
has motivated new developments in instrumentation and 
simulations, the latter being the focus of this paper. 

Recent results obtained with amplitude (Michelson) 
interferometry have started to re v eal stars as extende d 
objects (e.g. iBaldwin etafl 1 19961 : iMonnier etafl l2007t ). 
and with non-uniform light intensity distributions in 
the milliarcsecond scale. Such interesting results can be 
further investigated with SII taking advantage of the 
longer (km) baselines and relative ease of observing at 
shorter (blue) wavelengths. For example, measuring stellar 
diameters at different wavelengths, will make it possible 
to further i nvestigate the wave length dependence of limb 
darkening, (|Mozurkewicbl I2003T ) and thus constrain stellar 
atmosphere models. Radii measurements with uncertainties 
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of a few percent, along with spectroscopic measurements 
are necessary to c onstrain the p osition of stars in the 
HR diagram (e.g. iKervellal l200Sl ). With the methods 
described in this paper, we show that diameters can in 
principle be measured with accuracies better than 1% when 
using realistic array configurations for future experiments 
such as CTA (Cherenkov Telescope Array). As another 
example we can consider fast rotating B stars, which are 
ideal candidates for i maging ob l atene ss, pole brightening 
l|Monnier et alj 20071: von Zeipej [l9 24t) . radiatively driven 
mass loss IIFriend fc Abbottl I 198b ). and perhaps even 



pulsation modes 1 Saio et all 120061 ) . The impact of rotation 
on stellar evolution is non-triv ial, and several studies hav e 
been made in the subject (e.g. lMartinlll996l ; lMaededll997T ). 
Images of rotat ing stars have become available in the pas t 
few years (e.g. IMonnier et"aTll2007l : lAufdenberg et al.l l2006h 
and measurements of oblateness with accuracies of a few 
percent have been made. We will show that this is compa- 
rable to what can be achieved with SII using large arrays of 
Cherenkov telescopes. There is also the case of interacting 
binaries, for which we can measure angular separation, 
diameters, and relative bright ness. It may even be possible 
to m easure mass transfer (jVerhoelst. van Aarle. fc Ackd 
2007). Measurements of the angular separation in binaries 
is crucial for determining the mas ses of stars. The se masses 
must be found to within ~ 2% (|Andersenlll99lf ) in order 
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to test main sequence models. With the methods described 
in this paper, we show that angular separations can be 
found to within a few percent from reconstructed images. 
A more convenient and accurate analysis (beyond the scope 
of this paper) does not require imaging, and should yield 
sub-percent uncertainties, so that testing main sequence 
models will become possible with SII. 



In preparation for a large-scale SII observatory de- 
ployment, several lab oratory experiments are in progress 
l|LeBohec et~aH l20ld ). Their main goal is to measure 
light intensity correlation between two receivers. It is also 
worth mentioning the Star Base observatory (located in 
Grantsville, Utah) which consists of two 3m light receivers 
separated by 24 m and which will be used to test high time 
resolution digital correlators, and to measure the second 
order degree of coherence for a few stars. Various analog 
and digital correlator technologies |Dravins et al.l ll994T ) 
are being implemented, and cross correlation of streams of 
photons with nanosecond-scale resolution has already been 
achieved. 



Intensity interferometry, unlike amplitude interferome- 
try, relies on the correlation between intensity fluctuations 
averaged over the spectral band at electronic (nanosecond) 
time resolution. These averaged fluctuations are much slower 
than the (femtosecond) light wave period. This correlation 
is directly related to the complex de gree of coherence 7^ as 
l|Labevrie. Lipson fc Nisensonll2006l ) 



< AZ, Alj > 

< hXlj > ' 



(1) 



Here, < h > is the time average of the intensity 
received at a particular telescope i, and Ah is the intensity 
fluctuation. Measuring a second-order effect results in 
lower signal-to- noise ratio when compa red to amplitude 
interferometry l|Holder fc LeBohed 120061 ). This sensitivity 
issue can be dealt with by using large light collection areas 
(such as those available with air Cherenkov telescopes), 
longer exposure times and baseline redundancy. The low 
frequency fluctuation can be interpreted classically as the 
beat formed by neighboring Fourier components. Since SII 
relies on low frequency fluctuations, which are typically 
several orders of magnitude smaller than the frequency 
of optical light, it does not rely on the phase difference 
between light waves, but rather in the difference between 
the relative phases of the two components at the detectors 
l|Hanburv Brownll 19741 ). The main advantage is the relative 
insensitivity to atmospheric turbulence and the absence of 
a requirement for sub-wavelength precisi on in the optics 
and delay lines (Han burv Brown] 1 1974! ). The complex 
mutual degree of coherence 7 is proportional to the Fourier 
transform of the radiance distribution of the object in the 
sky (Van Cittert-Zernike theorem). However, since with 
SII, the squared-modulus of 7 is the measurable quantity, 
the main disadvantage is that the phase of the Fourier 
transform is lost in the measurement process. The loss 
of phase information poses a severe difficulty, and images 
have in the past been reconstructed from the bispectrum 
techn ique, using monolithic apertures (e.g. lLawrence et al.l 
1990). The imaging limitations can be overcome using a 
model-independent phase recovery technique. Even though 



several phase reconstruction techniques exist (|Fienudll98lT ). 
we concentrate on a two dimensional version of the one 
dimen sional analysis introduced by iHolmes fc Belen'kiil 
(2004), which is based on the theory of analytic functions, 
and in particular the Cauchy-Riemann equations. 

Following recent successes in Gamma ray astron- 
omy, a next generation Cherenkov telescope array is in a 
preparatory stage. This proje ct is currently known as CTA 
(Cherenkov telescope array) ()CTA Consortium! 12010). and 
will contain between 50 and 100 telescopes with apertures 
ranging between 5 m and 25 m. In this paper we investigate 
the sensitivity and imaging capabilities of SII implemented 
on such an atmospheric Cherenkov telescope array. We start 
with a discussion of sensitivity (section [2}, followed by a 
discussion of simulating noisy data as would be realistically 
obtained with such an array (section [3} . Since data have a 
finite sampling in the (u, v) plane, we discuss our method of 
fitting an analytic function to the data in order to estimate 
derivatives which are needed for phase reconstruction (sec- 
tion [4]). We then briefly discuss phase recovery (section [5]l 
and proceed to quantify the reconstruction quality using sev- 
eral criteria. We start with the simple case of uniform disks 
(section 16. ip and progressively increase the degree of im- 
age complexity by including oblateness (section l6.2|l . binary 
stars (section 16. 3[) . and obscuring disks and spots (section 
El. 



2 SENSITIVITY 

The signal to noise ratio (SNR) for an intensity correlation 
measurement depends on the degree of correlation 7, the 
area A of each of the light receivers, the spectral density 
n (number of photons per unit area per unit time, per fre- 
quency), the quantum efficiency a, the electronic bandwidth 
Af j and the observation time t. The SNR can be expressed 
as (jHanburv Brownll 19741 ) 



SNR = n(A, T, m v ) A a 7 2 ^ Aft/2. 



(2) 



The spectral density n is related to the visual mag- 
nitude m, of the star as well as its temperature T and 
observing wavelength A. The spectral density n(A, T, m v ) 
is the number of black body photons per unit area, per 
unit frequency and per unit time. The dependence of the 
visual magnitude m v is found by recalling that the flux 
for a th magnitude star with a temperature of 9550° K 
obser ved at 550 nm is 3.64 x 10" 23 Wm" 2 Hz" 1 l|Besselll 
Il979t ). This in turn corresponds to a spectral density of 
10~ 4 m~ 2 s _1 Hz _1 . The spectral density as a function of 
temperature (for different visual magnitudes and observing 
wavelengths) is shown in Fig. [T] and we see that at con- 
stant visual magnitude and observing wavelength, higher 
temperatures correspond to higher spectral densities. We 
find that the increase in temperature AT(A,T, Am v ) is 
approximately xk ^ Am v for the range of temperatures and 
wavelengths considered in Fig. [T] For example, at 400 nm, 
a decrease of 1 visual magnitude is equivalent to increasing 
the temperature of the star from 5000 K to 5700 K. For this 
reason high temperature objects are easier targets for SII. 
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Figure 1. Spectral density as a function of temperature for sev- 
eral different visual magnitudes and observed wavelengths. Atmo- 
spheric absorption as a function of the wavelength is not taken 
into account. 
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We use a preliminary design of th e CTA project as 
an array configuration ( Be rnloml l2008h . which is shown 
in Fig. [2] This array contains N — 97 telescopes and 
N(N — l)/2 = 4646 baselines (many of which are redun- 
dant) which are shown in Fig. [3] Each detector is assumed 
to have a light collecting area of 100 m 2 and a light detection 
quantum efficiency of a = 0.3. Using a X/D criterion, we 
find that the largest baselines of 1.5 km resolve angular 
scales of ~ 0.05 mas at 400 nm. The smallest 48 baselines 
of 35 m resolve angular scales of ~ 2 mas . However, we 
show in section 16.11 that the largest angular scales that 
can be realistically imaged with our analysis, in a model 
independent way, are more determined by baselines of 
~ 70 m. This is because the estimation of derivatives of the 
phase (needed for phase recovery) degrades as the number 
of baselines is reduced. Baselines of 70 m resolve angular 
scales of ~ 1.2 mas at 400 nm. 

These order-of-magnitude considerations are taken into 
account when performing simulations and image reconstruc- 
tions, i.e. the minimum and maximum size of pristine im- 
ages that can be reconstructed by data analysis, do not go 
far beyond these limits. More precise array resolution lim- 
its are presented in section 16.11 (diameters ranging between 
0.06 mas — 1.2 mas) . By combining these angular scales 
with the SNR (eq. we obtain Fig. U This figure shows 
the highest visual magnitude, for which photon correlations 
(with |7| = 0.5) can be detected (5 standard deviations), 
as a function of the temperature, and for several different 
exposure times. Also shown in Fig. [4] is the shaded region 
corresponding to angular diameters between 0.03 mas and 
0.6 mas Q and observable within 100 hrs. From the Figure 

1 These curves of constant angular size can be found approx- 
imately by recalling that the visual magnitude m v is related 
to a calibrator star of visual magnitude mo by: (m v — mo) = 
—2.5 log F/Fq. Here, F and Fo refer to the flux in the visual band 
. To express m v — mo as a function of the angular size, note that 
flux is proportional to (9 2 T 4 , where 9 is the angular size and T is 
the temperature of the star. 



Figure 2. Array configuration used for our analysis. 
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Figure 3. A total of 4656 non-zero baselines are available in the 
array design used in this study. Gray scale measures the log of 
the baseline redundancy. Since the array shown in Fig.[2]is almost 
symmetric with respect to x and y, only a quadrant of the (u, v) 
plane is displayed. 

we can see how correlations of photons from faint stars can 
be more easily detected if they are hot. To quantify the num- 
ber of stars for which photon correlations can be detected 
with the IACT a rray, we perform use the JMMC stellar di- 
ameters catalog (Lafras se et al ] l20ld ). We find that ~ 1000 
(out of ~ 33000) stars from the JMMC catalog can be de- 
tected within 1 hr, correlations from ~ 2500 stars can be 
detected within 10 hrs, and ~ 4300 can be detected within 
50 hrs. In Fig. [2] we show a random sample of 2000 stars 
(out of ~ 33000) from the JMMC catalog. Interstellar red- 
dening may play a role in reducing the number of measurable 
targets 
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Figure 4. The four parallel curves indicate the maximum de- 
tectable visual magnitude as a function of the temperature for 
several exposure times. Each of these four curves corresponds to 5 
standard deviation measurements and \^\ = 0.5. Also shown is the 
(shaded) region corresponding to angular radii between 0.03 mas 
and 0.6 mas, and observable within 100 hrs of observation time. 
The positions in (T, m v ) space of 2000 star s from the JMMC 
stellar diameters catalog ( Lafr asse et alj|201Cl ) are included. 

3 SIMULATION OF REALISTIC DATA 

The original "pristine" image consists of 2048 x 2048 pixels 
corresponding to ~ 10 mas x 10 mas of angular extension 
and a wavelength of A = 400 nm. The Fourier transform 
of the image is then performed via an FFT algorithm and 
normalized so that its value is one at zero baseline. This 
results in a Fourier transform sampled every ~ (8 m) /A, 
i.e. 2 x 10 7 cycles per radian field-of-view at a wavelength 
A of 400 nm. We then find the squared-modulus of the 
degree of coherence between the members of each pair of 
telescopes. This is obtained from a linear interpolation 
of the squared Fourier magnitude in the finely sampled 
FFT. Diurnal motion is not taken into account in the 
simulations. Diurnal motion plays a significant role in 
increasing the [u, v) coverage when exposure times are 
long. As a consequence there is less (u, v) coverage in the 
simulations since projected baselines do not drift with time. 
The smaller (u, v) coverage is however compensated by 
smaller statistical error in the correlation measurements. 

The final step in the simulation phase is the addition 
of noise to the correlation at each baseline. This noise was 
found to be Gaussian by performing the time integrated 
product of two random streams of simulated photons as 
detected by a pair of photo-multiplier tubes. The standard 
deviation of the noise added to each pair of telescopes is 
calculated from eq. [2] In this paper we take the signal 
bandwidth to be A/ = 200 MHz. An example of simulated 
data as a function of telescope separation is shown in Fig. 
[5] This corresponds to a 3 rd ma gnitude uniform disk star 
(T = 6000°K) of radius 0.1 mas and 10 hrs of observation 
time. The software used for the simulations, as well as the 
analysifl, was developed in C. 



2 See sections|4]and[5]for details on the analysis. Some variants of 
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Figure 5. Example of simulated data for a 3 rd ma gnitude uni- 
form disk star (T = 6000° K) of radius 0.1 mas and 10 hrs of ob- 
servation time. Here we show 7 (instead of the directly measured 
|7| 2 ) as a function of telescope separation. 



4 FITTING AN ANALYTIC FUNCTION TO 
THE DATA 

The estimation of derivatives of the Fourier log-magnitude 
is at the heart of the Cauchy-Riemann phase recovery 
algorithm (section [SJ), and is thus greatly simplified when 
data is known on a square grid rather than in a 'randomly' 
sampled way as is directly available from observations. 
Once simulated data are available (or observations in the 
future), an analytic function is fitted to the data. 

An analytic function can be expressed as a linear com- 
bination of basis functions. When data f(xi) = |7(a;i)| 2 are 
available at baselines Xi, with uncertainty Sf(xi), the coeffi- 
cients of the basis functions can be found by minimizing the 
following x 2 '- 

2 _ (f( x i) ~ J2k akgk{aR(xi))) ] 

x -2-^ Sf(x t ) ■ (i) 

Each Ofe is the coefficient of a basis function g^- The 
constant a is a scaling factor, and R is a rotation operator. 
The scaling factor and rotation angle are found by first 
performing a two-dimensional Gaussian fit. Finding the 
appropriate scale and rotation angle has the advantage of 
reducing the number of basis elements needed to fit the data. 

Basis functions which tend to zero at very large 
baselines, where data are scarce (see Fig. El, are ideal. For 
this reason, we use Hermite functions. There are situations 
where data are more easily fit with a different set of basis 
functions, e.g. a binary with unresolved members. In such 
a situation, data do not rapidly tend to zero at large 
baselines, so the Hermite function fit may contain a large 
number of elements and result in high frequency noise where 



the analysis software were developed in MATLAB. All software 
is available upon request. 
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data are scared- The best choice of basis functions may 
therefore depend on the structure of the object. However, 
for consistency, we use the Hermite fit for all the objects 
that we analyze, and find that it gives reasonably good 
results. 

The x 2 minimization problem can be turned into a lin- 
ear system by setting the set of partial derivatives to 
zero. We typically start with a small number of basis el- 
ements, say eight, and only increase the number of basis 
elements if the optimized reduced \ 2 ls greater than some 
predefined value. 



5 CAUCHY-RIEMANN PHASE 
RECONSTRUCTION 

In order to perform a model-independent image re- 
construction, the phase of the Fourier transform needs 
to be recovered from m a gnitud e information only 
l|Labevrie. Lipson fc Nisensonl l200ut ). Since the image 
is real, the Fourier magnitude is symmetric with respect 
to the origin in the [u, v) plane. Lack of phase informa- 
tion results in the reconstructed image being arbitrary 
up to a translation and reflection. We discuss the one- 
dimensional phase reconstruction prior to developing the 
two-dimensional analysis. We can start by first approxi- 
mating the continuous Fourier transform I(x) by a discrete 
one, i.e. I(mAx) = 0{jAe)e iimk ° AxAe , where O{0) 
is the image in the sky and ko is the usual wave vector. 
Then it becomes convenient to express the discrete Fourier 
transform in phasor representation , i.e. I(z) = R(z)e t *^ 
where z = e tmk o AxA0 j s complex. Since the discrete Fourier 
transform is a polynomial in z, the theory of analytic 
functions can be applied. As a consequence, we obtain the 
Cauchy-Riemann equations in polar coordinates, which 
relate the log-magnitude and the phase along the purely 
real and imaginary axis. If data were available along the 
purely real or purely imaginary axis, we could directly 
solve for the phase by integrating along these directions. 
However, since z, the independent variable of the Fourier 
transform (z = e *™*aA*A«^ hag modulus equal to 1; the 

phase differences that we seek lie along the unit circle in the 
complex plane. On the other hand, one can show that by 
using the Cauchy-Riemann equations, the phase differences 
along the radial direction in the complex plan^fl are directly 
related to the differences in the logarithm of the magnitude 
which is available from the data (see appendix A for more 
details) . 

The procedure to find the phase consists of first 
assuming a general solution form, then taking differences 
in the radial direction of the complex z-plane, and finally 
fitting the data to the radial differences of the assumed 
solution. A general form of the phase can be postulated by 
noting that the phase is a solution of the Laplace equation 

3 A basis set consisting of products of sines and cosines is more 
suitable in this situation 

4 If ^ is the real axis and ip is the imaginary axis, then a difference 
along the radial direction is A£ + iAi/j (assuming A£ = a£ and 
Aip = aiji, where a is a proportionality constant). 



in the complex plane (applying the Laplacian operator 
on the phase and using the Cauchy-Riemann equations 
yields zero). Since the phase differences are known along 
the radial direction in the complex plane we can take 
radial differences of the general solution and then fit the 
log-magnitude differences (available from the data) to the 
radial differences of the general solution. 

We can think of this one-dimensional reconstruction as 
a phase estimation along a single slice in the Fourier plane. 
A generalization to two dimensions can be made by doing 
the same procedure for several slices as described in Fig. [5] 
In fact, the requirement that a two-dimensional complex 
function (z x ,z y ) be analytic, is equivalent to s atisfying the 
Cauc hy-Riemann equations in both z x and z y (|Hormanderl 
1990). The direction of the slices is arbitrary, however for 
simplicity we reconstruct the phase along an arbitrary set 
of perpendicular directions in the Fourier plane, and noting 
that one can relate all slices through a single orthogonal 
slice, i.e. once the phase at the origin is set to zero, the 
single orthogonal slice sets the initial values for the rest of 
the slices. 

One can also require that the phase at a particu- 
lar point in the complex plane be exactly equal when 
reconstructed along z x or z y since each reconstruction 
is arbitrary up to a constant (piston) and a linear term 
(tip/tilt). However, imposing this requirement results in a 
severely over-determined linear system. More precisely, by 
imposing equality in n 2 points in the complex plane, and 
having 2n slices (each with an unknown constant and linear 
term) , results in a linear system of n 2 equations and 4n un- 
knowns. Alternative methods of requiring slice consistency 
are a possible way of improving phase reconstruction, but 
are beyond the scope of this paper. 

The Cauchy-Riemman approach, with horizontal or 
vertical slices, and a single orthogonal slice, gives reason- 
ably good results, however, it is not the only possible ap- 
proach. We have also investigated Gerchberg-Saxton phase 
retrieval, Generalized Expectation Maximization, and other 
variants of the Cauchy-Riemann approach. It is prema- 
ture to conclude which of these approaches is best at this 
time, given the limited imagery and SNR levels that have 
been explored. However, the Cauchy-Riemman approach 
has shown to give better resul ts in a number of cases 
|Holmes. Nunez fc LeBohedl2010f ). 



6 IMAGING CAPABILITIES 

We investigated the imaging capabilities for simple object^], 
namely uniform disk-like stars, oblate rotating stars, bina- 
ries, and more complex images. For most of the objects that 
we consider, image reconstruction is not necessary, i.e. from 
the Fourier magnitude alone, one can extract radii, oblate- 
ness, relative brightness in binaries, etc. Estimation of these 
parameters is probably more accurate when extracted di- 
rectly from Fourier magnitude data only, especially if some 

For preliminary study see N unez et al. 
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Figure 6. Schematic representation of two-dimensional phase re- 
construction approach. Several parallel slices are related to a sin- 
gle orthogonal slice. 



a-priori knowledge of the original image is available. How- 
ever, measuring simple parameters from reconstructed im- 
ages is the first step in quantifying reconstruction capabil- 
ities with I ACT arrays. We assume no a-priori knowledge 
of the images that are being reconstructed, and then do a 
statistical study of the uncertainties of the reconstructed 
parameters. 



6.1 Uniform disks 

In order to quantify the uncertainty in the reconstructed 
radius, we simulate data corresponding to 6 th magnitude 
stars (T = 6000° K) with disk radii up to 1 mas for 50 hours 
of exposure timfl An example of such a reconstruction is 
shown in Fig. [7] where the radiance is shown in arbitrary 
units between and 1. For a uniform disk, the reconstructed 
phase is null in the first lobe, and we find that the rms 
deviations from the true phase are approximately 0.19 rad 
in the null zone. A first look at the reconstruction in Fig. [7] 
reveals that the edge of the reconstructed disk is not sharp, 
so a threshold in the radiance was applied for measuring 
the radius. The radius is measured by counting pixels 
above a threshold and noting that the area of the disk is 
proportional to the number of pixels passing the threshold. 
After experimenting with different radii, we chose the 
threshold for measuring the radius to be 0.2. We can now 
compare the simulated and reconstructed radii as is shown 
in Fig. |8j where each point in the Fig. corresponds to an 
individual simulation (including noise) and reconstruction. 
Further optimization in the threshold for measuring the 
radius should yield a slope even closer to unity in Fig. [8] 

Fig. [8] clearly shows that stellar radii ranging from 
0.03 mas to 0.6 mas can be measured with uncertainties 
ranging between sub-percent and a few percent (Figures 
and I10p . It can be seen from Fig. [8] that the uncertainty 



6 This brightness and exposure time correspond to uncertainties 
in the simulated data of a few percent. Such long exposure times 
can be reduced to a few hours as is shown in Fig. 1101 



Figure 7. Example of a reconstructed uniform disk of radius 
0.1 mas. Also shown is a slice of the reconstructed image (solid 
line) compared to a slice of the pristine image convolved with the 
PSF of the array (dashed line). 
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Figure 8. Simulated vs. Reconstructed radii for magnitude 6 
stars with 50 hours of observation time (see footnote O . The top 
sub-figure shows the uncertainty for a 0.2 mas measurement. The 
bottom sub-figure shows the residual (Reconstructed-Real) along 
with the uncertainty in the radius. 



increases linearly as a function of the pristine (simulated) 
radius. This is due to a decrease in the number of baselines 
that measure a high degree of correlation when the angu- 
lar diameter increases. As the pristine radius 9 decreases, 
the distance to the first zero in the correlation increases as 
~ 9 , so the number of telescopes contained within the 
airy disk increases as ~ 9~ . Consequently, decreasing the 
pristine radius is equivalent to increasing the number of in- 
dependent measurements by a factor of ~ 9~ 2 . Since the 
uncertainty decreases as the square root of the number of 
independent measurements, the error decreases linearly with 
the radius. For radii above 0.6 mas, there are simply not 
enough baselines to constrain the Fourier plane information 
for image reconstruction. For radii greater than 0.6 mas, the 
distance to the first zero in the degree of correlation is of 
the order of 100 m, but only baselines at 35 m and 50 m are 
capable of measuring the Fourier magnitude with more than 
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Figure 9. Histogram of real radius minus reconstructed radius for 
50 hours of exposure time on a 6 th magnitude star (T = 6000 K) 
of 0.1 mas radius. 
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Figure 11. Simulated and reconstructed oblate rotator of mag- 
nitude 6 and 50 hours of observation time. 
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Figure 10. Percent error as a function of time for two recon- 
structed radii. This error was found by performing several recon- 
structions for each radius and exposure time, and then taking the 
standard deviation of the reconstructed radius. 



Figure 12. Real vs. reconstructed a/b for oblate stars. The dis- 
tance between the two linear envelopes is 2 standard deviations. 
All pristine images that have either a > 0.5mas or b > 0.4mas 
are not included 



3 standard deviations (see eq.[2}. In Fig.[l0]we show the rel- 
ative percent error (RMS of a statistic) as a function of time 
for two different radii, where it can be seen that a relative 
error of a few percent is achieved after only a few hours. 

6.2 Oblate stars 

For oblate stars we use the same magnitude and exposure 
parameters that are used for disk-like stars. Uniform oblate 
stars can be described by three parameters: the semi-major 
axis a, the semi-minor axis b, and the inclination angle 0. 
Judging from the limitations obtained from reconstructing 
disks, we produce pristine images whose values for a 
and b are random numbers less than 1 mas, and choose 
1 ^ a/b ^ 2 . The value of the inclination angle 9 also 
varies randomly between 0° and 90°. A typical image 
reconstruction can be seen in Fig. 1111 

After applying a threshold on pixel values as was done 
for disk shaped stars, the reconstructed parameters are 
found by calculating the inertia tensor of the reconstructed 
image. The eigenvalues and eigenvectors of the inertia 
tensor provide information for the reconstructed values 



of a, b and 6. To do this, we make use of the relation 
between the matrix eigenvalue and the semi- major/minor 
axes I X x = \a 2 M, where M is the integrated brightness. A 
similar relation for I yy holds for the semi-minor axis b. 

The resulting reconstructed semi-major/minor axes 
as a function of their real values have a similar structure 
as the scatter plot for reconstructed radii shown in Fig. 
[8] and are well reconstructed up to 0.5 mas within a few 
percent. In Fig. 1121 it can be seen how the uncertainty in 
the reconstructed oblateness a/b increases with increasing 
oblateness. As with disk shaped stars (section 16. ip . the 
uncertainty in the reconstructed semi- major/minor axes 
decreases as the square root of the number of baselines 
measuring a high degree of correlation. Therefore, the 
uncertainty in the reconstructed semi- major /minor axes is 
proportional to ~ \fab~, and the error in the reconstructed 
oblateness is proportional to ~ y^a/b + a 3 /b 3 . 

The reconstructed inclination angle as a function of the 
pristine angle is shown in Fig. 1131 and several factors play 
a role in the uncertainty of the reconstructed value. For a 
fixed value of a and b, the orientation of the telescope ar- 
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Figure 13. Real vs. reconstructed inclination angle for oblate 
stars. All pristine images that have either a > 0.5 mas, b > 0.4 mas 
or a/b < 1.1 are not included. We also note that reconstructed 
angles are always less than 45° due to the fact that the Fourier 
magnitude data is the same for pristine images flipped about the 
x, y or x and y axes. Therefore, for all pristine angles larger than 
45°, we replace the reconstructed angle by S' rec = 90° — Orec- 



ray with respect to the main lobe of the Fourier magnitude 
determines the number of baselines that measure a high de- 
gree of correlation. The number of baselines that measure a 
high degree of correlation is greater when the main lobe of 
the Fourier magnitude is aligned with the x or y direction 
of the array (see Fig. [2]), and is smaller by a factor of ~ \/2 
(assuming a uniform grid of telescopes) when its main axis 
is at 45° with respect to the array. However, the uncertainty 
(proportional to spread of points) in Fig. [13] does not ap- 
pear to be symmetric at 0° and 90°, and is smaller at 90°. 
This due to the way the phase is reconstructed, i.e. due to 
the slicing of the Fourier plane along the u or v directions 
(see section [5J. In the case of Fig. 1131 the (u, v) plane is 
sliced along the u direction, with a single orthogonal refer- 
ence slice along the v direction. The main lobe of the Fourier 
magnitude of an oblate star has more slices passing through 
it when it is elongated along the v direction (corresponding 
to an inclination angle of 90° in image space), yielding a 
better reconstruction. This is in contrast to the orthogonal 
case of 0°, where the main lobe of the Fourier magnitude 
has a smaller number of slices passing through it. 



6.3 Binary stars 

Simulated data are generated for 5 th magnitude binary 
stars (T = 6000° K), and an exposure of 50 hours after 
noting that the uncertainty in the degree of correlation is 
of the order of a few percent (eq. Binaries stars are 
parameterized by the radii ri and V2 of each star, their 
separation d, position angle 9, and relative brightness in 
arbitrary units between and 1. We generate pristine 
images with random parameters within the following 
ranges: radii are less than 0.3 mas, angular separations 
are less than 1.5 mas, the relative brightness per unit 
area is less than or equal to 1, and the orientation angle is 
less than 90°. A typical reconstruction can be seen in Fig. 1141 
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Figure 14. Simulated and reconstructed binary of magnitude 6 
and 50 hours of observation time. 
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Figure 15. Real vs. reconstructed angular separation in binary 
stars. Binary stars whose relative brightness is less than 0.3 are 
not included in this plot. 



old of 0.2. We then find the radius for each bright spot and 
its centroid position. Identifying spots is a non trivial task 
in noisy reconstructed images and our analysis sometimes 
fails to identify the "correct" reconstructed spots. For 
example, a common issue is that close reconstructed spots 
that are faintly connected by artifacts, are sometimes 
interpreted as a single spot. It should be again stressed that 
image reconstruction may not be the best way to measure 
reconstructed parameters. For example, the data can just as 
well be fit by the general form of the Fourier magnitude of 
a resolved binary system (containing only a few parameters) . 

In Fig. 1151 we show reconstructed angular separations 
as a function of their real values. The reconstructed values 
of the angular separation are found to within ~ 5% of their 
real values and cannot be much larger than what is allowed 
by the smallest baselines. We find that stars separated by 
more than d m ax ~ 0.75 mas are not well reconstructed since 
the variations in the Fourier magnitude start to become 
comparable to the shortest baseline. 



To measure the reconstructed parameters we identify 
the two brightest spots whose pixel values exceed a thresh- 



In Fig. [16] we show the reconstructed values of the 
radii as a function of their pristine values. We find ~ 10% 
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Figure 16. Real versus reconstructed radii for binary stars. We 
only include cases where one of the members is less than 3 times 
brighter than the other. 



uncertainties in each of the reconstructed radii. Aside 
from the angular separation, a variable that plays a role 
in successfully reconstructing pristine radii is the ratio of 
absolute brightnesfl of both binary members. When one 
of the two members is more than ~ 3 times brighter than 
the other, the fainter star is found to be smaller than its 
pristine value, and sometimes not found at all when one of 
the members is more than ~ 10 times brighter than the 
other. This is in part because the sinusoidal variations in 
the Fourier magnitude start to become comparable to the 
uncertainty. For example: a non-resolved binary star with 
one component 20 times brighter that the other, has relative 
variations of ~ 10%. With all the redundant baselines, a few 
percent uncertainty in the measured degree of correlation is 
sufficient to accurately measure these variations. However, 
when the binary components are resolved, the relative 
variations decrease with increasing baseline and baseline 
redundancy is not sufficient to reduce the uncertainty in the 
measurement of the Fourier magnitude. This signal to noise 
issue can of course be improved by increasing exposure time. 

There are also issues related to algorithm performance. 
One such problem has to do with the fit of the data to an 
analytic function (see section [4|. When the scale of the fit 
(found by an initial Gauss fit) is found to be too small, too 
many basis elements are used to reconstruct the data, and 
high frequency artifacts appear in reconstructions. Small 
initial scales are typically related to the binary separation 
as opposed to the size of individual components, and it is 
the latter which correctly sets the scale of the fit. Artifacts 
may be then mistaken for binary components, and incorrect 
reconstructed parameters may be found. Results improve 
significantly when either the correct scale is set or (model- 
independent) image post-processing is performed (see sec- 
tion [63]). 



6.4 Featured images 

We now present two examples of more complex images, and 
show that the capabilities can, to a large extent, be under- 
stood from results of less complex images, such as uniform 
disks and binaries. In Fig. [17] we show the reconstruction of 
a star with a dark band (obscuring disk), corresponding to 
a 4 th magnitude star and 10 hrs of observation time. The 
metric used to quantify the agreement with the pristine 
image (bottom left corner of Fig. I17[) is a normalized 
correlatiorjfl whose absolute value ranges between (no 
correlation) and 1 (perfect correlation/anti-correlation). To 
quantify the uncertainty in the correlation, we perform the 
noisy simulation and reconstruction several times, and find 
the standard deviation of the degree of correlation. In the 
case of Fig. [17] the correlation c is c = 0.947 ± 0.001. Note 
that the uncertainty in the correlation is apparently small, 
and the image reconstruction is not perfect, which implies 
that the reconstruction is not only affected by the SNR 
level, but also by the reconstruction algorithm performance. 

In order to determine the confidence with which 
we can detect the feature (dark region within the disk), 
this correlation is compared to the correlation of the 
reconstruction, and a featureless image, i.e. a uniform disk 
whose radius matches the radius of the pristine image. This 
comparison allows us to quantify the confidence with which 
we actually detect the feature (dark spot). Consequently, 
the correlation of the reconstructed image with a uniform 
disk is c = 0.880 ± 0.001, which is lower by 61 standard 
deviations. 

Another example of a complex image reconstruction 
can be seen in Fig. 1181 Here we show the reconstruction 
of a star with a dark spot, whose correlation with the 
pristine image is c = 0.940 ± 0.001. We can compare this 
correlation with the correlation of the reconstructed image 
and a uniform disk, which is c = 0.904 ± 0.001 (lower by 30 
standard deviations). 

For both examples, we also calculate the correlation 
with the pristine image as a function of the angular scale 
(in mas) of the pristine image. We then find the degree of 
correlation of each reconstruction and its pristine image, 
and also the degree of correlation of the reconstruction with 
a uniform disk. By comparing these two correlation values, 
we can find the smallest feature that can be reconstructed. 
Below some point we no longer distinguish between the 
reconstruction of the featured image and that of a uniform 
disk. We find that the smallest feature that can be recon- 
structed is close to 0.05 mas. This can already be seen from 
the order of magnitude estimate made in section [2] and a 
comparable value of 0.03 mas is found in section IHTl When 
pristine images have angular sizes greater than 0.8 mas, the 
degree of correlation drops significantly due to lack of short 
baselines. Note that this value is higher than what is found 
in section |6~T1 This is due to higher signal to noise ratio in 



8 For two images A% and Bij, the normalized correlation Cj 

7 Product of relative brightness per unit area (in arbitrary linear ls — Max^^^N (<ta&b) ^i,j^i+k,j+l 1 where 

units between and 1) and relative area of both stars a a and <jb are the standard deviations of images A and B. 
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Figure 17. Star with obscuring disk (raw reconstruction). This 
corresponds to 4 th magnitude and 10 hrs of observation time. 
The correlation between the real and reconstructed image is c = 
0.947 ± 0.001. Note that an inverted gray scale is used. 



Reconstructed Ima 
1 




Figure 18. Star with dark spot (raw reconstruction). This cor- 
responds to a 4 th magnitude star, 10 hrs of observation time and 
a degree of correlation of 0.940 ± 0.001. 
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Figure 19. Simulated and (raw) reconstructed oblate rotator of 
magnitude 3 and 10 hours of observation time. Reconstruction 
quality is lower when compared to Fig. 1111 This is due to high 
frequency noise in the Hermite fit of the data, and also due to 
limitations of the phase reconstruction algorithm. Post-processing 
is applied and Fig. 1201 is obtained. 
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Figure 20. Image post-processing on Fig. ll9l The MiRA software 
was used to optimize on an area 0.3 mas X 0.3 mas. 



the simulated data corresponding to Figures IPTl and [T8l The 
resolution limits discussed in turn correspond to ~ 16 x 16 
effective resolution elements (pixels). 



6.5 Post-processing 

Image post-processing is performed in order to improve the 
raw images presented above. The type of post-processing 
that is currently being investigated is analogous to the 
data analysis performed in conventional amplitude interfer- 
ometry. In this analysis, the image is slightly modified so 
as to maximize the agreement between the data and the 
magnitude of the Fou r ier tr ansform of the reconstructed 
image. See iThiebauti (|2009l ) for details. Also, one can 
introduce additional and very general constraints to this 
maximization procedure, but this is beyond the scope of this 
paper. Such an analysis depends strongly on the starting 
image, which can be provided by the analysis presented in 
this paper. An example of image post-processing performed 



on the oblate star of Fig. [19] is shown in Fig. 1201 

We have performed this type of post-processing on 
the images in F igures [T71 and 1181 using the MiRA software 
(|Thiebautl l2009), and we obtain the ones shown in Figures 
1211 and [221 Aside from optimizing agreement with the data, 
no additional constraints are imposed. These preliminary re- 
sults show the overall reduction in noise and improvement 
in the sharpness of the reconstructed images. A systematic 
study of the improvements with image post-p rocessing is 
currently being investigated (|Nunez et al.ll20Tlh . 



7 CONCLUSIONS 

We perform a simulation study of the imaging capabilities 
at 400 nm of an IACT array consisting of 97 telescopes 
separated up to a km. This is a preliminary design for 
the CTA project, expected to be operational in 2018. Our 
method uses a model-independent algorithm to recover 
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Figure 21. Post processed Fig.[T3 The MiRA l|Thiebaudl200gh 
software was used to optimize on an area of 1 mas X 1 mas. 
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Figure 22. Post processed Fig. 1 181 The MiRA software was used 
to optimize on an area 0.3 mas X 0.3 mas. 



the phase from intensity interferometric data. We test the 
method on images of increasing degrees of complexity, 
parameterizing the pristine image, and comparing the 
reconstructed parameter with the pristine parameter. 
We now summarize our results and briefly discuss how 
fundamental stellar parameters can be constrained with the 
methods described above. 

We first find that for bright disk-like stars (M v — 6), 
radii are well reconstructed from 0.03 mas to 0.6 mas. 
Even though using a Cauchy-Riemann based approach 
to recover images might not be the most efficient way 
to measure stellar radii, such a study starts to quantify 
the abilities of measuring other scale parameters in more 
complicated images (oblateness, distance between binary 
components, etc.). The range of angular radii that can 
be measured with a CTA-like array (0.03 — 0.6 mas) 
will complem ent existing measurements (2 — 50 mas) 
[|Hanifj|200lh . With the aid of photometry, the effective 



temperatur^fl scale of stars within 0.05 — 0.5 mas can 
be extended. Multi-wavelength angular diameter mea- 
surements will als o reveal the wavelength dependenc e 
of limb-darkening (|Aufdenberg. Ridgwav fe White! 120101 ). 
and such measurements can be used to constrain energy 
transport models as is done in ampl itude interferometry 
l|Aufdenberg. Ludwig fe Kervellall20 



For oblate stars, similar results for extracting geomet- 
rical parameters are found. Due to the relative ease of SII 
to observe at short (~ 400 nm) wavelengths, measuring 
fundamental parameters of hot B type stars is possible. 
B stars are particularly interesting since rapid rotation, 
oblateness, and mass loss are a common feature. We show 
that oblateness can be accurately measured, and a next 
step will be to quantify the capabilities of imaging realistic 
surface brightness distributions in hot stars. By imaging 
brightness (temperature) distributions, we will be able 
to study effects s uch as limb darkenin g and mass loss in 
hot massive stars ilRidgwav et alj|2009h . as well as gravity 



not massive stars (ttidgwav 
darkening (Ivon Zeipellll924l ) 



Binary stars are well reconstructed when one of the 
members is not much brighter (three times as bright) 
than the other, and when they are not too far apart 
(^ 0.75 mas). As with amplitude interferometry, SII, along 
with spectroscopy, will allow the determination masses and 
orbital parameters i n binary stars. I f measured with enough 
precision « 2%) (| Andersen! Il99lh , the determination of 
the mass can be used to test main sequence stellar models. 
An advantage of using an array such as the one used in 
this study, is that individual radii can be resolved. An 
interesting phenomena to b e studied with int eracting binary 
stars is mass transfer (e.g. IZhao et al.ll2008h , and capabili- 
ties for imaging this phenomena can be further investigated. 

Two examples of reconstructions of more complex im- 
ages are presented, and demonstrate the capability to resolve 
features in the sub-milliarcsecond scale. Results show im- 
provement when post-processing is performed. An optimiza- 
tion procedure improves the image reconstructions found by 
the methods presented above by maximizing their agreement 
with the data. A more comprehensive study of reconstruc- 
tions of complex images is currently being performed. 



APPENDIX A. 
RECOVERY 



CAUCHY-RIEMANN PHASE 



Recall that since SII can only allow to measure \j\ 2 , the 
phase information is lost, and we would like to recover it for 
imaging purposes, with only the magnitude information. If 
we denote I(z) = R(z)e ti '^ z \ where z = £+iip, we obtain the 
following relations from the Cauchy-Riemann equationif^l: 



Defined as T, 



L 

y 4irR 2 c 



,1/4 



where R is the radius, L is the 



luminosity, and a is the Stephan-Boltzmann constant. 
10 The C-R equations can be applied because "I" is a polynomial 
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9$ 

dip 



dlnR = 
dlnR 

d~f 



ds 



ds_ 



(4) 
(5) 



where we have defined s as the log-magnitude. Notice 
the relation between the magnitude and the phase. By us- 
ing the Cauchy-Riemann equations we can write the log- 
magnitude differences along the real and imaginary axes as: 



ds . j. 9$ . . 
A ^ = d^=^ 
ds 9$ 



(6) 

(7) 
(8) 



If the log-magnitude were available along purely the £ 
or the ip axes, we could solve the previous two equations for 
the phase. 

However, notice that because \z\ — 1, we can only 
measure the log-magnitude on the unit circle in the complex 
space (£, ip). 

In general, we can write the log-magnitude differences 
along the unit circle as 



Note from eq. (|10[1 that the length in the complex plane 
associated with Asm is A<p = | A£+iAip\, and that the length 
associated with A$x is Ap — \A^+iAtp\ . Now setting p — 1, 
Ap — A<f>, and for simplicity of presentation, expanding for 
small A(j>, we obtain 



A4>x(<£) = j A(f> (a.j cosj(p + bj sin jq 



(14) 



equations (|9I10[1 from the measured Ash, and thus 3> can be 
found in the complex plane, with an uncertainty in ao and bo. 
The coefficients a,j and bj can be calculated by performing 
the following integrals: 

Oj = - — r / — —cos j <j>d(j) 



2nj 
1 



(15) 



(16) 



Note however that the previous expressions must exist, 
which is not the general case. More explicitly, if the magni- 
tude is zero, the log-magnitude is singular. When imaging 
finite objects in image space, there will always be zeros in 
the magnitude of the Fourier transform. In practice we are 
always sample limited and nothing prevents us from calcu- 
lating cij and bj approximately. 



ds . _ ds . , 
d~^ + ^ 
9® . „ 9$ . , 

A$L ■ 



(9) 
(10) 



Here A$x corresponds to phase differences along a 
direction perpendicular to Asm, that is, perpendicular to 
the unit circle in the £ — ip plane. We are however interested 
in obtaining A3>ii, so that we can integrate along the unit 
circle. 

The general form of $ can be found by taking second 
derivatives in eq. ((5| and thus noting that $ is a solution of 
the Laplace equation in the complex plane. 



<9 2 $ _ 
Si 5 " + ^ ~ 



(11) 



The general solution of <f>(z) in polar coordinates (p,4>) 
is (|jacksonlll998h 

$(p, cj>) = a + b 4> + p J (oj cos j<j> + bj sin jtj>) , (12) 

3 

where terms singular at the origin (p~ J ) have been 
omitted. These singular terms lead to ambiguous recon- 
structions including flipped images and have not been found 
to be essential for most reconstructions. 

Now taking the difference of $ along the radial direction 
we obtain 

A$ x (p, 0) = V pP((l+^Py-l) ( aj cos j4> + bj sin j<t>) .(13) 
is p 



REFERENCES 

Andersen J., 1991, AAPR, 3, 91 

Aufdenberg J. P., Ludwig H. G., Kervella P., 2005, ApJ, 
633, 424 

Aufdenberg J. P. et al., 2006, ApJ, 645, 664 

Aufdenberg J. P., Ridgway S., White R., 2009, astro2010: 

The Astronomy and Astrophysics Decadal Survey, 2010, 

8 

Baldwin J. E. et al, 1996, A&A, 306, L13 
Bernlohr K., 2008, American Institute of Physics Confer- 
ence Series, 1085, 874 
Bessell M. S., 1979, PASP, 91, 589 
CTA Consortium, 2010, 1 arXiv: 1008.37031 
Dravins D. et al., 1994, Proc. SPIE 2198, 289 
Fienup J. R., 1981, JOSA, 71, 1641 
Friend D. B., Abbott D. C, 1986, ApJ, 311, 701 
Hanbury Brown R., 1974, The Intensity Interferometer, 

Taylor & Francis, London. 
Haniff C. A., 2001, Galaxies and their Constituents at the 

Highest Angular Resolutions, 205, 288 
Holder J., LeBohec S., 2006, ApJ 649, 399 
Holmes R. B., Belen'kii M. S., 2004, JOSA 21, 697 
Holmes R. B., Nunez P. D., LeBohec S., 2010, Proc. SPIE 
7818 

Hormander L., 1990, An introduction to complex analysis 
in several variables. Elsevier Science 

Jackson J. D., 1998, Classical Electrodynamics, 3rd Edi- 
tion, by John David Jackson, Wiley- VCH. 

Kervella P., 2008, IAU Symposium, 252, 405 

Labeyrie A., Lipson S. C, Nisenson P., 2006, An Intro- 
duction to Optical Stellar Interferometry, by A. Labeyrie, 
S. G. Lipson and P. Nisenson, Cambridge University 
Press, Cambridge 



High angular resolution imaging with stellar intensity inter •ferometry using air Cherenkov telescope arrays 13 



Lafrasse S., Mella G., Bonneau D., Duvert G., Delfosse X., 

Chelli A., 2010, VizieR Online Data Catalog, 2300, 
Lawrence T. W., Fitch J. P., Goodman D. M., Massie N. A., 

Sherwood R. J., 1990, Proc. SPIE, 1237, 522 
LeBohec S. et al., 2010, Proc. SPIE, 7734, 7734 
Maeder A., 1997, A fcA, 321, 134 
Maier G., preprint (|arXiv: 0907.5 1 18) 
Martin E. L., Claret A., 1996, A&A, 306, 408 
Monnier J. D., 2007, Bulletin of the American Astronomi- 
cal Society, 38, 854 
Monnier J. D. et al., 2007, Science, 317, 342 
Mozurkewich D. et al, 2003, ApJ 126, 2502-2520 
Nunez P. D., Lebohec S., Kieda D., Holmes R., Jensen H., 

Dravins, D., 2010, Proc. SPIE 7734 
Nunez P. D. et al., 2011, In prep. 

Ridgway S. et al., 2009, astro2010: The Astronomy and 

Astrophysics Decadal Survey, 2010, 247 
Saio H. et al, 2006, ApJ, 650, 1111 
Thiebaut E., 2009, New Astron. Rev. 53, 312 
Verhoelst T., van Aarle E., Acke B., 2007, A&A, 470, L21 
von Zeipel H., 1924, MNRAS, 84, 665 
Zhao M. et al., 2008, ApJ, 684, L95 



